Heavy Quark Diffusion with Relativistic Langevin Dynamics 
in the Quark-Gluon Fluid 
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The relativistic diffusion process of heavy quarks is formulated on the basis of the relativistic 
Langevin equation in Ito discretization scheme. The drag force inside the quark-gluon plasma 
(QGP) is parametrized according to the formula for the strongly coupled plasma obtained by the 
AdS /CFT correspondence. The diffusion dynamics of charm and bottom quarks in QGP is described 
by combining the Langevin simulation under the background matter described by the relativistic 
hydrodynamics. Theoretical calculations of the nuclear modification factor _Raa and the elliptic 
flow V2 for the single electrons from the charm and bottom decays are compared with the experi- 
mental data from the relativistic heavy ion collisions. The Raa for electrons with large transverse 
momentum (px > 3 GeV) indicates that the drag force from the QGP is as strong as the AdS/CFT 
prediction. 

PACS numbers; 24.85.+p, 05.40.Jc, 11.25.Tq 



I. INTRODUCTION 

The physics of the quark-gluon plasma (QGP) is ac- 
tively studied by means of the relativistic heavy-ion 
collisions at Relativistic Heavy Ion Collider (RHIC) in 
BNL and will be pursued further at Large Hadron Col- 
lider (LHC) in CERN The space-time evolution of 
the heavy-ion collisions at RHIC is well described by 
the (3-|-l)-dimensional relativistic hydrodynamics sup- 
plemented with the hadronic cascade after chemical 
freczcout Q. Information on the collective dynamics of 
QGP is obtained by the soft probes such as distributions 
of Hght hadrons at low momentum, while the information 
of microscopic dynamics of QGP is obtained by the hard 
probes such as jets, heavy quarks, and heavy quarkoni- 
ums 0. 

In the present paper, we focus on charm and bottom 
quarks which behave as impurities in QGP. Experimen- 
tally, the signal of the heavy quarks can be extracted from 
the single electron spectra through semileptonic decays 
Theoretically, the energy loss of heavy quarks in 
QGP has been estimated in perturbative QCD (pQCD) 
techniques [1, 0|- However, it was pointed out recently 
that the convergence of the weak coupling expansion of 
the drag force for heavy quarks is rather poor at the 
temperature relevant to RHIC and LHC, so that the cal- 
culation in the leading order would not be reliable Q. 
Possible alternative way to estimate the drag force is to 
use the duality conjecture between the gauge theory and 
string theory (AdS/CFT correspondence) Al- 
though it can be applied only to the A/" = 4 supersymmet- 
ric Yang-Mills plasma with large 't Hooft coupling, the 
result obtained may provide us with a hint for the drag 
force in the strong coupling regime of the QCD plasma 
if appropriate translation is made [l^ . 

The purpose of this paper is to study the connection 
between the drag force acting on the charm and bot- 
tom quarks in QGP and the final electron spectra. To 
make such connection, we introduce relativistic Langevin 



equation for heavy quarks under the background of the 
quark-gluon fluid described by the ideal hydrodynamics. 
We need relativity since the transverse momentum of the 
heavy quarks at RHIC is not necessarily smaller than 
their rest masses. Our relativistic Langevin equation is 
formulated in Ito discretization scheme. The diffusion 
constant and the drag force are related through a gener- 
alized fluctuation-dissipation relation. As for the drag 
force, two distinct models, pQCD and AdS/CFT, are 
considered. To calculate the space-time dynamics of light 
quarks and gluons, (3-|-l)-dimensional hydrodynamics is 
used, which is necessary to calculate the electron spectra 
of different impact parameters in the heavy ion collisions. 
The Langevin equation for heavy quarks is numerically 
solved from the initial distribution generated by Monte 
Carlo generator PYTHIA ^ until the freczcout of the 
heavy quarks. The transverse- momentum (px) depen- 
dence of the nuclear modiflcation factor (Raa) and the 
eUiptic flow (W2) of single electrons as decay products 
of heavy quarks are calculated and compared with the 
RHIC data. 

This paper is organized as follows. In Sec. [Hi after for- 
mulating the relativistic Langevin equation and a gener- 
alized fluctuation-dissipation relation, we introduce two 
extreme models of the drag force motivated by pQCD 
and AdS/CFT. In Scc. lIIIl the relativistic hydrodynamics 
for light degrees of freedom and the relativistic Langevin 
equation for heavy degrees of freedom are combined in 
order to describe the heavy quark diffusion in dynam- 
ical QGP fluid. The initial condition of heavy quarks, 
the algorithm of Langevin simulation in dynamical back- 
ground, and the treatment of the freezeout and decay of 
heavy quarks are discussed in detail. In Sec. IIVI the nu- 
merical results of our calculation are presented. We show 
the proflle of heavy quark diffusion, heavy quark spectra, 
and single electron spectra and compare our single elec- 
tron spectra with experimental data. Physical meaning 
of our results are also discussed. Section V is devoted 
to summary and concluding remarks. In Appendix A, 
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we show a derivation of the relativistic Kramers equa- 
tion from the relativistic Langevin equation in the Ito 
discretization scheme. 



II. LANGEVIN DYNAMICS OF HEAVY 
QUARKS 

In this section, we formulate the relativistic Langevin 
equation in the local rest frame of the fluid. A gener- 
alized form of the fluctuation-dissipation relation is de- 
rived. Then, we discuss the drag forces calculated in the 
pQCD approach and in the AdS/CFT approach. Finally, 
wc introduce a phcnonicnological model of the drag co- 
efficient r and the diffusion constant D which satisfy the 
generalized fluctuation-dissipation relation. 



A. Relativistic Brownian motion 

Suppose that there exist a time scale tb during which 
a Brownian particle changes its momentum and a micro- 
scopic time scale Tm during which light dynamical de- 
grees of freedom change state and lose time correlation. 
If these time scales satisfy tb ^ Tm, one can describe the 
Brownian particle by the Langevin equation [l^ . In such 
a situation that the charm/bottom diffuses inside the 
quark-gluon plasma, there are some extra complications: 
(i) The background quark-gluon fluid expands rapidly in 
space and time with the local 4- velocity u^{x, t), and (ii) 
the initial momentum distribution of the charm/bottom 
governed by the hard QCD process has high momentum 
component larger than their quark masses. 

As for (i), we define the Langevin equation in the rest 
frame of matter {u^ = (1,0,0,0)) and the real motion 
of the Brownian particle is obtained by the local Lorentz 
boost back to the moving frame. As for (ii), we take 
into account the relativistic kinematics of the Brownian 
particle in a minimal way by using relativistic dispersion 
relation E{p) ~ x/p^ + . Then the Langevin equation 
in the rest frame of matter with minimum relativistic 
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kinematics may be written as [15 1 



-At 



Am = -r{p)pAt+at). 



(1) 

(2) 



Here Ax{t) = x(t') - ^{t), and Ap{t) = v{t') - P{t), and 
At = t' — t \s Si discrete step of time. The momentum 
dependent drag coefficient is denoted by r(p) which is re- 
lated to the time scale of the Brownian particle tb ~ F"^. 
Also ^{t) is a noise obeying the probability distribution 
VF[^(t)] which we take to be the Gaussian white noise 
with a normalization constant C: 



W[£,{t)]d?£.{t) = Ccxp 



2D{p)At 



d'm- (3) 
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(4) 
(5) 



where D{p) is a momentum dependent diffusion constant. 
Note that ^(t) is not a single microscopic kick but a sum 
of microscopic kicks during the time At. 

Throughout this paper, we use the Ito discretization 
scheme of the Langevin equation; namely, all the argu- 
ment of p in the right hand side of Eqs. ^ and ^ are 
evaluated at the pre-point time t. This is particularly 
useful for numerical simulations due to obvious reason. 
The relativistic Kramers equation, which is a partial dif- 
ferential equation for the probability of the particle dis- 
tribution in the phase space P{p,x,t), is then obtained 
as (see the derivation in Appendix A) 



d_ 

dp 



T{p)p 



ld_ 

2W 



D{p))p{p,x,t). 



(6) 



Note that 8/ dp acts not only on P but also on F(p) and 
D{p). 

Demanding that Eq. ([6|) has the relativistic Maxwell- 
Boltzmann distribution (the Jiittner distribution) 
P{p, X, t) (X exp[— -y/p^ _|_ M'^/T] as a stationary solution, 
wc obtain a constraint between the drag and the diffusion 
as 



T{p) + G{p) = 



D{P) 
2ET' 



(7) 



with G{p) = dD{p)/2pdp = dD{p)/d{p^). If D is p- 
indcpendcnt, Eq. ([7]) reduces to the relativistic analogue 
of the Einstein relation F = tMtp = ^ „P.„ obtained in 

r L Zrjl rj Zl\l 1 



B. Modeling the energy loss of heavy quarks 

Energy loss of heavy quarks in the deconfined phase 
has two sources; the collisional energy loss due to elas- 
tic scattering of a heavy quark with the plasma con- 
stituents and the radiative energy loss associated with the 
induced emission of the gluon. In the leading order (LO) 
of the weak-coupling QCD perturbation, these processes 
are found to have different momentum dependence of the 
heavy quark and could become comparable in magnitude 
[1, 0] . Recently, the convergence of such weak-coupling 
expansion was questioned by an explicit calculation of the 
collisional process in the next-to- leading order (NLO) 
The drag coefficient for the 3-color, 3-flavor QCD in the 
non-rclativistic kinematics (il/ ^ T, p) reads 



F 



pQCDU/»T,p 



M 



(- In <7-H 0.07428 -FL 



(8) 
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For the QCD coupling constant relevant at RHIC and 
LHC {g ~ 2), the weak-coupling expansion has an obvi- 
ous problem of convergence. From the phenomenological 
point of view, it has been argued in the past that rel- 
atively large drag force is necessary to account for the 
RHIC data 

Alternative approach to the drag force is provided by 
the AdS/CFT correspondence In the A/" = 4 

super- Yang-Mills theory (SYM) at large Nc and large 't 
Hooft coupling A = g'^-^^N^, energy loss of an external 
quark with velocity v is obtained as 



dp 
'dt 



7r\/A 2 



2 "sY^Tf 



rp-2 f 

^ SYM ~ ■ 



M 



(9) 
(10) 



Here the first equation is valid for arbitrary mass of ex- 
ternal quark, while the second equation is valid for M 3> 
\f\T [Tl| . By matching the energy density and the heavy 
quark potential in the SYM plasma to those in QCD 
plasma, one finds Tsym — 7qcd/3^/^ and 3.5 < A < 8.0 
|13 |. Then, the drag coefficients may be estimated as 



TAds/CFT (2.1 ± 0-5)^- 



(11) 



A remarkable feature of this formula in contrast to the 
weak-coupling estimate is that F is p- independent [ll|. 
The fundamental question here is, of course, the reliabil- 
ity of the translation from SYM to QCD both conceptu- 
ally and numerically. 

Given the theoretical uncertainties in estimating the 
drag force as mentioned above, we will take a phenomeno- 
logical approach in this paper: We adopt the paramet- 
ric dependence of the drag coefficient motivated by the 
AdS/CFT in Eq. (HIl) with the overaU magnitude left as 
a free parameter: 



(12) 



The dimensionless drag coefficient 7 is assumed to be in- 
dependent of r, M, and -p throughout this paper. The 
corresponding diffusion constant D is obtained from the 
generalized fluctuation-dissipation relation in Eq. ([7]) 
with the physical boundary condition, D ^ as F ^ 0: 



D = 2ET • F • 1 



2T^ 

= 1^{E + T). (13) 



It is in order here to make two remarks on the dy- 
namics we employed in Eqs. p2)) and (fT3l) . (i) Since 
we have assumed F to be p-independcnt motivated by 
AdS/CFT, the diffusion constant D depends necessarily 
on the momentum of the Brownian particle. One may 
alternatively assume that D is independent of p while F 
depends on p as F(p) = D/{2E{p)T) [ll|. Such dynam- 
ics would simulate the p-dependence of the drag force due 
to coUisional process in the weak-coupling regime Q . (ii) 



At ultra-high energies p ^ M , the dominant energy loss 
occurs through the induced emission of the gluons. In 
this case, the condition tb ^ is violated. Thus the 
Langevin approach becomes inapplicable [ll[ and a dif- 
ferent approach based on radiative energ y lo ss is required 
to describe heavy quarks in the QGP 0, UMl ■ With these 
reservations in mind, we consider our ansatz Eqs. (|12p 
and (|13p as phenomenological but characteristic dynam- 
ics of QCD and try to estimate 7 from the observed single 
electron data at RHIC in later sections. 



III. HYDRO + HEAVY-QUARK MODEL 
A. Background QGP fluid 

The hydrodynamics has been quite successful in the de- 
scription of collective flow phenomena in heavy ion colli- 
sions at RHIC. Since the hydrodynamics gives space-time 
evolution of temperature and flow velocity of the fluid so 
that the local rest frame of fluid is well-defined. Then, 
the Langevin equation in the previous section formulated 
in the local rest frame of the fluid is applicable directly. 

Let us first summarize the relativistic hydrodynamic 
model ^, tlB, i2(i, i21,] whose basic equation reads 



df.T'"' = 0. 



(14) 



Here T'^'^ is energy-momentum tensor. For strongly in- 
teracting matter with zero viscosity, T^'' becomes 



T"" = {e + P)u>'u'' - Pg' 



til' 



(15) 



where e, P, and u'^ arc energy density, pressure, and four 
fiuid velocity, respectively. The baryon chemical poten- 
tial is neglected, because it is small near mid-rapidity 
at RHIC energies. We solve Eq. in the Bjorken 

and 



coordinates (r, x, y, 77^), where t 



ijs = 5 lii[(t -|- -j)/(t — 2)] are proper time and space-time 
rapidity, respectively. 

In the high temperature (T > = 170 MeV) QGP 
phase, we employ the bag equation of state (EOS) for 
massless partons (u, d, s quarks and gluons) with i?^/* = 
247.19 MeV. Here the bag constant is tuned to have tran- 
sition to the hadron resonance gas at Tc- In the hadron 
phase {T < Tc = 170 MeV), a resonance gas of hadrons 
with the mass up to A(1232) is employed [2l|. Volume 
fraction of QGP /qgp in the mixed phase is 



/qgp 



e - eiiad 
eqcp — ehad ' 



(16) 



where cqg? (ehad) is the maximum (minimum) value of 
the energy density in the mixed phase. Later we will 
utilize /qgp to define the effective lifetime of QGP and 
the freezeout condition for the heavy quarks. 

Hot QGP with local thermalization is assumed to be 
produced at tq = 0.6 fm. The entropy density distri- 
bution at To in the mid-rapidity is taken to be propor- 
tional to a linear combination of the number densities 
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of participants and binary collisions in the transverse 
plane [l^. For the initial condition of the flow veloc- 
ity, Bjorken's scaling solutioii, Ux(to) — Uy{To) =0 and 
Mz(to) = sinhTy^, is employed [22| . With these initial con- 
ditions, the hydrodynamic model can well reproduce the 
experimental data of charged particles at RHIC 

The space-time evolution of the QGP fluid obtained as 
above has been exploited for a quantitative study of hard 
and rare probes such as azimuthal jet anisotropy, nuclear 
modification factor of identified hadrons, disappearance 
of back-to-back jet correlation, J/tp suppression, and di- 
rect photon emission p3l |. 
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5.2 
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4.0 



TABLE I: (a) Relaxation times of charm and bottom quarks for 
7 = 0.3, 1.0, and 3.0 at T = 210 MeV. (b) Lifetimes of QGP for 
different centralities and freezeout conditions. In (a) Mc and Mb 
arc chosen to be 1.5 GeV and 4.8 GeV, respectively. In (b) we 
adopt two characteristic impact parameters b = 3.1 and 5.5 fm. 



B. Heavy quark diffusion in quark-gluon fluid 

We solve the Langevin equation ([T]) and ^ with the 
drag coefficient T given by Eq. p2)) in the local rest frame 
of the fluid element. The dimensionless parameter 7 is in- 
versely proportional to the relaxation time tq of a heavy 
quark as 



1 
f 



(17) 



The TQ is listed in Table IUja) for three typical values, 
7 = 0.3 (weak coupling), 7 = 1.0 (intermediate cou- 
pling), and 7 = 3.0 (strong coupling). The characteristic 
temperature felt by the heavy quark during the space- 
time history in the QGP fluid is taken to be 210 MeV (as 
for the reasoning of this number, see Sec. IIV AT|) . 

Let us now introduce an effective lifetime of QGP, 
TQGP; by the following definition: At r = tqgp, the QGP 
fraction /qgp in Eq. (|16[) a,tx^y~z = reaches to 
/o, which takes a value between and 1. For /o = 0, 
the effective lifetime is defined as the time when QGP 
disappears completely, while fo — ^ corresponds to the 
time when hadronic phase starts to appear. The effective 
lifetime of QGP is listed in Table |Jb) for two different 
impact parameters and for three different values of /q. 

From the comparison of tq and tqgp in Table HI one 
finds that the initial momentum distributions of charm 
and bottom quarks will be changed by QGP only slightly 
for the weak drag force (7 = 0.30). On the other hand, 
for the strong drag force (7 = 3.0), both charm and bot- 
tom quarks are affected by QGP and their momentum 
distributions would be modified substantially. 



1. Initial distribution of heavy quarks 

On the initial hypcrsurfacc tq = 0.6 fm, initial trans- 
verse positions of heavy quarks are distributed according 
to the overlap function of two nuclei A and B in the 



transverse plane Tab (2;, y): 

TABix,y) = Ta U+ Tb ) , 



TA{B){x,y) 



(18) 



where Pa(b) is the Woods-Saxon parametrization of nu- 
clear density. The heavy quarks are assumed to stream 
freely in the longitudinal direction for < r < tq and ac- 
quire the momentum rapidity yp ~ j]s- Thus the initial 
heavy quark distribution in the phase space reads 



dN 



(Pp(Px±TQdr]g 



d<^P?rj. I Arts 
-^TAB{x,y) 



To 



(19) 



The initial momentum spectrum of heavy quarks 
da^^ /d^p inp + p collisions is calculated by pcrturbative 
QCD to leading order (LO) utilizing the event generator 
PYTHIA 6.4 [il. 

In Fig. [1] the initial distribution of heavy quarks in 
the transverse plane and that in the momentum space 
at mid-rapidity {\yp\ < 1) are shown. The momentum 
distribution is normalized to the invariant cross section 
in p + p collisions. Note that the initial momentum dis- 
tribution of the charm quark has a steeper slope at high 
Pt than that of the bottom quark. Nuclear effects such 
as shadowing and Cronin effect are not considered for 
simplicity. 

In Fig. Ufa), we show the differential cross section of 
electrons from heavy quarks in p-\-p collisions obtained by 
PYTHIA (the LO perturbative QCD). Theoretical cross 
section underestimates the experimental value by a factor 
5-10 as shown in Fig. [DJb). The discrepancy is known 
to become smaller by taking into account higher orders 
beyond LO [l^, [2^ . Since we are mainly interested in the 
the "shape" of the px distribution above 3 GeV in this 
paper, we adopt the LO result for simplicity in spite of 
the problem of absolute magnitude in the LO calculation. 



2. Simulation of the Brownian motion 

The Langevin equation of a heavy quark is defined in 
the local rest frame of QGP. The information of local 
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(a) HQ configuration 
10 , ^ r 




3^ . f\-0 

::^i0 1 2345 6 789 10 
Pt [GeV] 



FIG. 1: (a) A sample of 500 heavy quarks at the initial in the trans- 
verse plane, for the Au + Au collision with impact parameter 5.5 
fm. (b) Invariant cross sections of charm and bottom production in 
p + p collisions in mid-rapidity {\yp\ < 1.0), which is proportional 
to the initial momentum distribution. 



flow velocity and local temperature at the position of 
the heavy quark is supplied from the relativistic hydro- 
dynamics. The algorithm of such Langevin simulation is 
summarized as follows. 

(i) Start from a sample of heavy quark at a position 
and a momentum according to the initial phase space 
distribution given in Eq. (|19p . 

(ii) Given the phase space location {p^,x^) in the 
laboratory frame, obtain the information of the local 
flow velocity u^{x) and local temperature T{x) from the 
output of the hydrodynamic simulation. 

(iii-a) Coordinate step: Make one discrete step for 
the heavy quark in the configuration space according 
to Eq. ([1]) by using discrete proper-time step As = 
{M/E)M: 



(a) c+b->e* in p+p 



• Data 
-LO pQCD 




20 

D 

^ 15 

Q. 

2 10 h 



< 

Q 



'(b) Ratio ^ 



10 



Pt [GeV] 



FIG. 2: (a) Experimental cross section for electron production in 
p + p collision at mid-rapidity [23 and the leading order pQCD 
result by PYTHIA. (b) The ratio of the experimental data and the 
LO result. Theoretical calculations are performed at \yp\ < 0.35 
and then properly normalized to obtain the cross section. 



(iii-b) Momentum step: Move to the rest frame of the 
fluid element by the Lorentz transformation, p ^ k. 
Make one discrete step for the heavy quark in momentum 
according to Eq. ^ using As: 



Afc(s) 



-7 — TT^fcAs 



ij ^ss^ ' 



A/2 



a^), (21) 
E{E + T)As. (22) 



Aa;'^(s) 



M 



As. 



(20) 



Then, move back to the laboratory frame by inverse 
Lorentz transformation (fc + Afc — > p'). 

(iv) Repeat the steps (ii) and (iii) until the volume 
fraction of QGP in the mixed phase (/qgp) reaches /q. 

Several comments are in order here about this proce- 
dure. 

• We use the proper-time step As instead of the or- 
dinary time step At in our simulation, simply be- 
cause the former is a Lorentz scalar and thus easy 
to handle in going back and forth between the lab- 
oratory frame and the fluid rest frame. We choose 
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As = 0.01 fm in our simulation, which is much 
shorter than the relaxation time of the drag force 
parameter adopted in this paper. 

• Owing to the Ito discretization scheme, the momen- 
tum step in (iii-b) can be performed only by using 
the information of flow and temperature at the cur- 
rent position of the heavy quark in the phase space. 

• It is not clear whether we should stop the heavy 
quark diffusion at the point when the mixed phase 
starts to appear or at the point when the mixed 
phase disappears. We consider this uncertainty as 
a systematic error and consider the three cases as 
shown in Table [ijb), namely /o = 0, 0.5, and 1. 



3. Freezeout and decay 

Once the local temperature around the charm(bottom) 
quark becomes lower than Tc, it hadronizes into D (B) 
mesons. Since we need to calculate single electron spec- 
tra from the heavy quarks, we focus on the following 
semileptonic decays: Z? — > e for D decay, B —f e for 
primary B decay, and B D e for secondary B de- 
cay. The hadronization of heavy quarks and the decay 
of heavy mesons are calculated by using PYTHIA 6.4 
(l3l |. Since we employ independent fragmentation given 
by PYTHIA, the effect of quark recombination to form 
D or B mesons is not taken into account. Such simplifi- 
cation would be more reasonable for heavy quarks with 
higher transverse momentum. Therefore, it is the high 
Pt region (e.g. above 3 GeV) that is suitable to compare 
our results with the experimental data. 



C. Observables 

Medium modification factor Raa for single electrons is 
defined by 



-Raa(pt) 



1 dNA+A/dpT 

A^coii dNp+p/dpT '' 



(23) 



where A^coii is the number of binary collisions calculated 
from the Glauber model. Since the initial heavy quark 
distribution is assumed to be without nuclear effects and 
to scale as A^coU in our calculation, the deviation of Raa 
from unity is solely attributed to the heavy quark dif- 
fusion in the hot medium. The elliptic flow for single 
electrons is defined by 



V2{pt) 



I' 



' ''/V/ cos 20 



+ A 



(cos 20). (24) 



This quantity indicates how much momentum anisotropy 
around the collision axis is given to the heavy quarks from 
the background medium. 



IV. NUMERICAL RESULTS 

Before showing the numerical results in detail, let us 
first summarize the basic parameters of our simulation. 
(1) The dimensionless drag coefficient 7 is a parameter 
to control the diffusion of heavy quarks in QGP. We take 
three characteristic values, 7 = 0.3, 1.0, and 3.0 cor- 
responding to weak, intermediate, and strong coupling, 
respectively. (2) The impact parameter b controls the 
volume and the lifetime of QGP. Thus it affects indi- 
rectly the heavy quark spectra at their freezeout and the 
single electron spectra. In all of the figures below except 
for Fig. ^ b is taken to be 5.5 fm (10-20% centrality). 
(3) The criterion of stopping the heavy quark diffusion 
in the mixed phase is given by /o which takes a value 
between and 1. Its dependence on the final results is 
considered to be a systematic error of our calculation. In 
all of the figures except for Fig. [9l we show the results at 
the central value, fo = 0.5. 



A. Heavy quark spectra 

1. Profile of heavy quark diffusion 

To estimate how long a heavy quark stays in the QGP 
region in terms of local fluid proper time, we define the 
"stay time" as 

ts = J2 = H As(^/A/)|frf 

steps steps 

= ^ As(p.VA/)|LF, (25) 



steps 

where FRF and LF imply the fluid rest frame and labo- 
ratory frame, respectively. By averaging over the heavy 
quarks starting initially with p™ and ending in mid- 
rapidity (|j/p| < 1.0) at their freezeout, we obtain the 
average stay time (<s)- 

Shown in Fig. [3] is the averaged stay time of heavy 
quarks as a function of their initial transverse momen- 
tum. The diffusion coefficient is taken to be 7 = 0, 0.3, 
1.0, 3.0, and 30.0. Here 7 = corresponds to the free 
streaming. On the other hand, 7 — 30.0 corresponds to 
the extremely strong coupling where the relaxation times 
at typical temperature 210 MeV are 0.22 fm for charm 
and 0.72 fm for bottom: The initial information on px is 
completely lost after a few fm of diffusion in this case. 

The figure shows that, for heavy quarks with large 
initial velocity compared to the fluid velocity (p^'" > 
1 GeV, p^'" > 3 GeV), the stay time becomes shorter 
for higher px because they get out of the medium in 
shorter times. Also, as the drag force becomes stronger, 
the stay time becomes longer as expected. As for the 
heavy quarks with small initial velocity (p^"^ < 1 GeV, 
p^'" < 3 GeV), the stronger the drag force, the shorter 
the stay time, since the drag force from the background 
fluid accelerates them more strongly. 
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FIG. 3: (Color online) The averaged stay time (tg) of (a) charm 
quarks and (b) bottom quarks with the drag coefficient 7 = 0, 
0.3, 1.0, 3.0, and 30.0 at mid-rapidity {\yp\ < 1.0). The impact 
parameter is chosen to be 5.5 fm in Au+Au collisions. For frcezcout 
condition, /o = 0.5 is adopted. 



FIG. 4: (Color online) The averaged temperature (T) of (a) charm 
quarks and (b) bottom quarks with the drag coefficient 7 = 0, 
0.3, 1.0, 3.0, and 30.0 in mid-rapidity {\yp\ < 1.0). The collision 
geometry and the freezeout condition are the same with those in 
Fig. |3] The fluctuation in high px in (a) is due to statistical errors 
of our simulation. 



Next we define the averaged temperature for the heavy 
quarks experienced during their stay in the QGP fluid: 

rEE(l/ts)^T(x)Ai|FRF. (26) 

steps 

By averaging over the heavy quarks starting initially with 
pJp and ending in mid-rapidity {\yp\ < 1.0) at frcezcout, 
we obtain the averaged temperature (T) shown in Fig. ID 

The figure shows that, for heavy quarks with large ini- 
tial velocity compared to the fluid velocity (p^p'" > 1 
GeV, p!;;'" > 3 GeV), the averaged temperature becomes 
higher for higher pip because they feel only the initial high 
temperature region before getting out of QGP. Also, as 
the drag force becomes stronger, the stay time becomes 
longer and averaged temperature becomes smaller. As 
for the heavy quarks with small initial velocity {p^™ < 
1 GeV, p^™ < 3 GeV), the stronger the drag force, the 
higher the averaged temperature, since they are strongly 
accelerated and quickly pass the low temperature region. 
It turns out that the average temperature lies between 
200 MeV and 220 MeV in the wide range of pip and 7; 
this is the reason why we adopted the typical tempera- 
ture 210 MeV in Sec. Ell 

Finally, let us define the the transverse momentum loss 



(momentum loss for short): 

Apr = pi^ - p°^^\ (27) 

where p^^ is the transverse momentum at the time of 
the freezeout of the heavy quark. By averaging over the 
heavy quarks starting initially with p]f and ending in 
mid-rapidity {\yp\ < 1.0) at freezeout, we obtain the av- 
eraged momentum loss (Apx) as shown in Fig. [5l 

For heavy quarks with larger initial pl^, the momen- 
tum loss per unit time (dynamical effect) is larger as 
seen in Eqs. ([2]) and ((T2)) while the average stay time 
(kinematical effect) is shorter. Therefore, there are two 
competing effects in the net momentum loss: In Fig. [SJ 
we find that larger initial momentum leads to larger mo- 
mentum loss, so that the dynamical effect wins over the 
kinematical effect. As for the dependence on drag coef- 
ficient, both dynamical and kinematical effects act ad- 
ditively for heavy quarks with large initial momentum 
(p^j^"^ > 1 GeV, pJj;'" > 1.5 GeV)) and the momentum 
loss is enhanced by increasing 7. For the heavy quarks 
with small initial velocity (p^j^™ < 1 GeV, p^J^'" < 1.5 
GcV), these two effects compete but we find in Fig. [5] 
that the dynamical effect seems to win, namely that the 
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FIG. 5: (Color online) The averaged momentum loss (Ap-p) of (a) 
charm quarks and (b) bottom quarks with drag coefBcicnts 7 = 
0.3, 1.0, 3.0, and 30.0 in mid-rapidity < 1.0). The collision 

geometry and the freezeout condition are the same with those in 
Fig. El 



FIG. 6; (Color online) R''^j^ of (a) charm quarks and (b) bottom 
quarks with drag coefficients 7 = 0.3, 1.0, 3.0, and 30.0 at mid- 
rapidity {\yp\ < 1.0). For collision geometry, we choose the impact 
parameter 5.5 fm in Au-(-Au collisions. For freezeout condition, the 
/o = 0.5 is adopted. 



stronger the drag force, the larger the momentum gain 
by the acceleration from the fluid. Note here that, for 
the extreme case 7 = 30.0, we have almost a linear in- 
crease of (Apx) = Pt ~ -Pt'* ^ function of pl^. This 
is simply because the heavy quarks are thermalized and 
p^p* is almost independent of p™. 



2. Nuclear modification factor -R^^ 

Let us define -R^^ {Q — c, b) for heavy quarks by re- 
placing the number of electrons Np^p (Na+a) in Eq- (ESI) 
by the number of heavy quarks at the freezeout Np_^_p 

(-^a+a)- This is a theoretical quantity not directly ac- 
cessible in experiment, but it is useful to examine the 
behavior of heavy quarks without the kinematical com- 
plication due to their semileptonic decays to electrons. 

Shown in Fig. [S] are R'aa f*-"^ charm and bottom in the 
mid-rapidity at impact parameter 5.5 fm as a function 
of p'^^ . There are two key factors which determine i? ! 
the momentum loss of heavy quarks and the initial distri- 
bution of heavy quarks. Starting from the initial distri- 
bution, the high momentum quarks loose energy due to 



drag force and are shifted to the low p?^*" region. There- 
fore, i?AA suppressed (enhanced) at high (low) p™*. 
This tendency is prominent for large drag force as ex- 
pected. Also, the suppression at high p™' is larger for 
the charm if we adopt the same 7. This is because the 
actual drag coefficient is ^T"^ /Mq so that the quark with 
smaller mass is affected more by the drag force. 



3. Elliptic flow «2 

In Fig. [71 we show the elliptic flow for the heavy quark 
V2 {Q ~ c, 6) at mid-rapidity at impact parameter 5.5 fm. 
It is clear that the charm and bottom quarks with any 
drag force at large p™* are less thermalized and thus they 
do not produce much momentum anisotropy. Note that 
the dominant contributions of heavy quarks with 7 = 
30.0 at large p™* may be those that start near the surface 
of the QGP fireball and with large initial p^f, therefore 
they are not much thermalized because of the too short 
stay times. On the other hand, charm quarks with small 
p'i^^ are thermalized for large drag force and develops 
V2 reflecting the flow of light particles. As for bottom 
quarks, they only have small momentum anisotropy with 
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FIG. 7: (Color online) of (a) charm quarks and (b) bottom 
quarks with drag coefficients 7 = 0.3, 1.0, 3.0, and 30.0 in mid- 
rapidity < 1.0). The collision geometry and the freezeout 
condition are the same with those in Fig. [6] In (a), the statistical 
errors for 7 = 3.0 and 30.0 are so largo at pJJ,"' > 6 GeV that we 
omit them. 



all drag forces but 7 = 30.0 at small because they 
are not enough thermalizcd. 




B. Electron spectra 
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1. Nuclear modification factor Raa 

Let us now examine the results of electrons and 
positrons (we call them just electrons for short) which 
arc the decay products from D and/or B mesons. In 
Fig. [5{a), (b), and (c), we show Raa of electrons (a) 
from charm quarks, (b) from bottom quarks, and (c) from 
charm-|-bottom quarks. The dependence of Raa on the 
drag coefficient 7 is understood easily: Larger drag co- 
efficient gives larger energy loss and Raa is suppressed. 
There is however one qualitative difference between Raa 
in Sec. IIV A 21 and Raa in the low pT region: Raa 
ceeds unity due to the shift of the high momentum quarks 
to low momentum quarks, while Raa stays around unity 
at low momentum. This is understood by recognizing 
that the low electrons come from wide range of heavy 
quarks with various freezeout momenta, so that low mo- 
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FIG. 8; (Color online) (a) Raa of electrons from charm, (b) -Raa 
of electrons from bottom, (c) -Raa of electrons from both charm 
and bottom, and (d) the ratio of electrons from the bottom and 
the net electrons. All results are in mid-rapidity (\yp\ < 0.35). The 
drag coefficient is taken to be 7 = 0.3, 1.0, and 3.0. The impact 
parameter is taken to be 5.5 fm in Au+Au collisions. For freezeout 
condition, the /o = 0.5 is adopted. In (d), the result of p -|- p 
collision calculated in the leading order pQCD by PYTHIA is also 
plotted. 
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FIG. 9; (Color online) Comparison of -Raa in our hydro + heavy- 
quark model with the experimental data 0, 13 • The Au-|- Au col- 
lision with the impact parameter (a) 3.1 fm and (b) 5.5 fm, both 
in mid-rapidity, \yp\ < 0.35. The drag coefficient is chosen to be 
7 = 0.3, 1.0, and 3.0 indicated by different colors. The freezeout 
condition is taken to be /o = 1.0, 0.5, and 0.0 which correspond 
to upper, middle, and lower points, respectively, within the same 
color. As for error bars in experimental data, we only plot the 
statistical errors 0il3. 



mentum electrons are not sensitive to the modification of 
the heavy quark spectrum due to diffusion. On the other 
hand, the electrons with high originate mainly from 
high heavy quarks and thus they are sensitive to the 
spectral change of heavy quarks. 

In Fig. \E[.d) , the number of electrons from bottom di- 
vided by that from charm-|-bottom for Au-|-Au collision 
is shown as a function of electron's pt together with that 
for p+p collision. In both p+p and A-l-A, more than 50% 
of electrons come from the bottom for px > 3 GeV. Fur- 
thermore, the ratio increases as the drag force becomes 
stronger. The kink structure of Raa at ^ 1 - 2 GeV 
in Fig. [5]^c) is understood by the fact that the dominant 
contribution to the electrons changes rapidly from the 
charm to the bottom. 

Finally we compare our numerical results with exper- 
imental data Q in Fig. [SI Here we show two cases of 
impact parameters 3.1 fm (0-10% centrality) and 5.5 fm 
(10-20% centrality) at mid-rapidity. The systematic er- 
rors due to the freezeout condition of heavy quark are 
represented by the three plots with the same color. Re- 



FIG. 10: (Color online) Comparison of 112 in our hydro -f- heavy- 
quark model with experimental data [3l in mid-rapidity {\yp\ < 
0.35). Experimental data of V2 is obtained in minimum bias anal- 
ysis, while our theoretical values of V2 are evaluated at impact pa- 
rameter 5.5 fm as a representative. The drag coefficient is chosen 
to be 7 = 0.3, 1.0, and 3.0 and the freezeout condition is /o = 0.5. 
As for error bars in experimental data, we only plot the statistical 
errors Q|. 



call that the comparison of our results and experimen- 
tal data is only reliable for pT > 3 GeV as discussed 
in Sec, nil B H and that bottom quarks are the dominant 
source of electrons in this region. 

Although definite conclusion cannot be made from the 
present comparison, it is likely that the intermediate to 
large value of the drag coefficient 7 = 1.0 - 3.0 is favored 
especially for small impact parameter. This number is 
rather close to the value 7 = 2.1 ± 0.5 predicted from the 
AdS/CFT correspondence (see Eq. (fTTj) ). We should re- 
mark, however, that the radiative energy loss 0, [13 and 
the relativistic diffusion via resonances combined with 
quark coalescence [l3| would be legitimate alternatives 
to describe the data, so that further systematic compar- 
ison of the data and theoretical calculations is called for. 



2. Elliptic flow V2 

We show our theoretical V2 of electrons in Fig. [10] as a 
function of px together with the experimental data 
Our i'2 does not depend much on the strength of the 
drag force for px > 3 GeV and stays small. Due to the 
poor statistic of both our simulation and the experimen- 
tal data in the relevant region, it is not clear whether 
theory and experiment are consistent with each other or 
not. Although it is still preliminary, recent PHENIX data 
show large V2 — 0.05 - 0.1 with small errors for 3 < pj < 5 
GeV at collisions with corresponding centrality [2y| . 



V. SUMMARY AND CONCLUDING REMARKS 

In this paper, we have examined the diffusion of 
heavy quarks in the dynamical QGP fluid on the ba- 
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sis of the rclativistic Langevin equation combined with 
the relativistic hydrodynamics. We establish a general- 
ized fluctuation-dissipation relation in Ito discretization 
scheme, Eq. ([7]), which relates the diffusion constant D{p) 
and the drag coefficient T[p) for the relativistic Brown- 
ian particle. Then we parametrized the drag coefficient 
motivated by the formula from the AdS/CFT correspon- 
dence for strongly coupled plasma, T = ^T'^ /M with the 
dimensionless coefScicnt 7 as a parameter. The space- 
time evolution of the QGP fluid composed of light quarks 
and gluons is treated by the full (3-|-l)-dimensional rcla- 
tivistic hydrodynamics for the perfect fluid. 

By solving the Langevin equation for heavy quarks 
under the influence of the QGP fluid, we obtain the 
space-time history of the diffusion process of charm 
and bottom in the realistic situation of the relativis- 
tic heavy-ion collisions. The initial space-time distribu- 
tions of charm/bottom are given by the event generator 
PYTHIA. The hadronization and the semileptonic decays 
of charm/bottom after they leave from the QGP region 
are treated by independent fragmentation and decay. Nu- 
clear effects for initial charm/bottom distributions and 
the quark recombination/coalescence in hadronization of 
heavy quarks, which would be important for low < 3 
GeV region of the flnal electron spectrum, are neglected 
for simplicity in this paper. 

Since we have the space-time history of the 
charm/bottom during their diffusion, we have looked at 
the average stay time of heavy quarks in QGP (ts), the 
average temperature felt by heavy quarks in QGP (T), 
and the average momentum loss (Ap) during the diffu- 
sion. We have also looked at the nuclear modification 
factor i?AA ^'^'^ elliptic flow of heavy quarks as a 
function of the transverse momentum of the heavy quarks 
at their freezeout p™*. The results indicate that, for suf- 
flcicntly large values of p™* > 3 GeV, there is a sizable 
suppression of i?AA large drag coefflcient, while one 
can see only a signiflcant effect in only for p™* < 3 
GeV which is not the region one can rely on our calcula- 
tion. 

Then we have compared our calculations of i?AA and 
the elliptic flow 1)2 for single electron with the RHIC data. 
First of all, the momentum distribution of the electrons 
do not necessary reflect the shape of the momentum dis- 
tribution of the heavy quarks at their freezeout due to 
decay kinematics. Also, the net electrons with > 3 
GeV are dominated by those from bottom quarks. A 
rough comparison of i?AA for px > 3 GeV suggests that 
the drag coefficient could be as large as 7 = 1.0 - 3.0. 
On the other hand, we are unable to extract useful in- 
formation from V2 for > 3 GcV because of the lack 
of statistics in both experiment and simulations. The 
value of 7 = 1.0 - 3.0 is consistent with that predicted 
by AdS/CFT approach for strongly interacting plasma 
(7 = 2.1 ± 0.5), although we could not exclude other 
descriptions of heavy quarks in QGP such as radiative 
stopping 0, HI] and the resonance scattering model [13] . 
High precision experimental data at RHIC and LHC for 



electrons from charm and bottom identified separately 
are highly desirable. Also the correlation of the trans- 
verse momenta of a heavy quark and a heavy anti-quark 
(and the associated electron-positron or D-D correlation 
[13, nil) could be a good observable to make detailed 
comparison of the theories and experiments. 

Before closing, we remark possible improvements of our 
approach to treat the region px < 3 GcV in a more reli- 
able way: (i) Initial heavy quark distributions beyond LO 
need to be considered for better control of their absolute 
magnitude, the px shape, and the charm/bottom ratio, 
(ii) nuclear effects on the initial charm/bottom distribu- 
tion need to be examined, and (iii) the hadronization 
of charm/bottom due to quark recombination processes 
needs to be taken into account. 
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APPENDIX A: 



RELATIVISTIC KRAMERS 
EQUATION 



In this Appendix, we derive the partial differential 
equation of P(p, x, t) or the Kramers equation using the 
Ito (prc-point) discretization scheme. We give here the 
general form of the relativistic Langevin equation as 



Af(t) 

{V^{tW)) 



Pit) 



E{p{t)) 
-A{p{i))p{i)At 

C ■ exp 



B{p{i))rf{t)At, 



At ' 



(Al) 



where E{p) ~ 



with Af being the mass of the 



Brownian particle. Here t = t corresponds to the Ito 
discretization and t = t + At corresponds to the Hanggi- 
Klimontovic discretization [l^. Also A{p), B{p), and 
fj{t) in the Ito discretization correspond to r(p), D(p), 
and i{t)/ D{p)At) in the text, respectively. Since the 
Langevin equation is based on Markovian process, one 
needs information only at time t' in order to know the 
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probability at later time t: 

P{p,x,t\po,xa,ta) (A2) 
d^p'd^x'P{p, X, t\jf, x', t')P{0, x\ i'lpo, 2-0, to), 



where t|Xo,io) {X = {p^x}) represents the condi- 

tional distribution function with a fixed initial condition 
Xq at time ip. In order to derive the partial differen- 
tial equation, we have to calculate P{X,t + At\X',t) 
from the Langevin equation. From the definition of 
P{X,t + At\X',t), 

P{X, t + M\X\ t) = {5{X - X{t + M))) \t,x' 
= {6[X ~ X' - AX{^{t),t)\) 

= ^([-AX(,Kt),i)]")^a^<5(X-X'), 

rn— 

(y(77(0,t)> ^ j d\{t)Wm]Y{r^{t),t). (A3) 

Here (• • • )\t,x' in the first line of Eq. (|A3|) represents 
the conditional probability with the fixed initial condition 
X{t) = X' . Note that in the last Hne of Eq. ((X3|l . the 
average is expressed by the variables at time t. Inserting 
Eq. ((X3|) into Eq. we obtain 

P(X,t + At|Xo,io) 
= y dX' P{X,t + At\X' ,t)P{X' ,t\XQM) 

J(X-X')+£([-AX(ry(t),t)]™) 

m— 1 

^^d^5{X - X')\- P{X\t\XM 
PiX,t\Xa,to) 

OO ^ 



dX' 



In the Ito discretization scheme, the relevant average 
values {[AX{rj{t),t)]"') are 



(Ax-it)) 

{Ap{t)) 
{Ap,{t)Ap^{t)) 



-At, 



m 

E{p{t))- 
-A{p{t))p{t)At, 
B{p{t))d,,At, 



(A5) 



and the others are in higher order in At. 

From Eqs. (jA4[) and (|A5p . the resulting relativistic 
Kramers equation reads 



d p d 
d 



P{p,x,t) 

1 d 
2 



(a(p)p + -—B{p))p{p, X, t). (A6) 



t\X„, h) + AtdtP{X, t|Xo, to)- 



(A4) 
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